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ABSTRACT 

Prompt 7-ray and early X-ray afterglow emission in gamma-ray bursts (GRBs) are characterized by a bursty 
behavior and are often interspersed with long quiescent times. There is compelling evidence that X-ray flares 
are linked to prompt 7-rays. However, the physical mechanism that leads to the complex temporal distribution 
of 7-ray pulses and X-ray flares is not understood. Here we show that the waiting time distribution (WTD) 
of pulses and flares exhibits a power-law tail extending over 4 decades with index ^ 2 and can be the man¬ 
ifestation of a common time-dependent Poisson process. This result is robust and is obtained on different 
catalogs. Surprisingly, GRBs with many (> 8) 7-ray pulses are very unlikely to be accompanied by X-ray 
flares after the end of the prompt emission (3.1 (t Gaussian confidence). These results are consistent with a 
simple interpretation; an hyperaccreting disk breaks up into one or a few groups of fragments, each of which 
is independently accreted with the same probability per unit time. Prompt 7-rays and late X-ray flares are 
nothing but different fragments being accreted at the beginning and at the end, respectively, following the very 
same stochastic process and likely the same mechanism. 

Subject headings: gamma-ray: bursts, waiting time distribution 


1. INTRODUCTION 

The first electromagnetic messenger of a gamma-ray burst 
(GRB) is the so-called 7-ray prompt emission, followed by 
the early X-ray afterglow on a timescale from minutes to 
hours. Long duration (> 2-3 s) GRBs are nowadays known to 
be associated with the core collapse of some kind of massive 
stars rid of hydrogen e nvelopes (see lWooslev & Bloomll200^ 
iHjorth & Bloomll2012l for reviews). Prompt 7-rays (with en¬ 
ergies in the keV-MeV range) are observed within a given 
GRB as a sequence of pulses (typically a few up to several 
dozens). In addition, for a sizable fraction of GRBs the fol¬ 
lowing decaying X-ray emission, which marks the end of 
the 7-rays, is characterized by the presence of X-ray flares 


which are sometimes observed as late as 10^ s (IBurrows et al. 

2005 

;IChincarini et al.l2007t 

Falcone et al.l2007HCurran et al. 

2008 

; iBernardini et al. 1201 1 

). Although mounting evidence 


exists that X-ray flares, like 7-ray pulses, result from the 
GRB inner engine activity rather than from external shocks 

(iLazzati & Pernal[2007l : IChincarini et al.ll2010l:lMargutti et af] 

I2OIOI) . key questions remain unanswered; what radiation pro- 

cess(es)? What information on the inner engine can we ex¬ 
tract? Is there a common process ruling inner engine activity 
across several decades in time? 

As a matter of fact, both emissions represent a tempo¬ 
ral point process, i.e. a time series characterized by the 
discrete occurrence of impulsive events superposed on a 
continuum. Intense bursting periods are often interspersed 
with relatively long (several up to tens of seconds) inter¬ 
vals with very low activity, compatible with the detector 
backg round, which are often referred to as quiescent times 
(QTs; i Ramirez-Ruiz & Merlo nill200lL Nak^ & PirarJ 120021 : 
lOuilligan et ani2002nDrago & Pagliaral2007l) . The study of 
the waiting time distribution (WTD), i.e. of how time inter¬ 
vals between adjacent peaks distribute, provides clues on the 
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nature of the stochastic process. 

In particular, it reveals the degree of memory and correla¬ 
tion and constrains the physical process responsible for the 
discontinuous and bursty release of energy. 

Processes showing similar on-off intermittency, or, equiv¬ 
alently, burs ty behavior or c lusterization, can be found in 
many fields (iPlatt et al.lll993h . The corresponding WTDs of¬ 
ten show power-law tails at long waiting times (WTs), whose 
index depends on the degree of clusterization of the time se¬ 
ries. Examples encompass the aftershock s equence ob served 
in earthquakes, described by Omori’s law (lUtsulllQ^ . neu¬ 
ronal firing activity, as well as a wide range of dynamical 
systems of human activity, such as mail and ema il exchanges 
( Oliveira & Barabasil2005tlEckmann et al.l2004l) . phone calls 
(iKarsai et al.ll2012 and reference s therein) all the way to vio¬ 
lent conflicts (iPicoli et al.ll2014ft . These processes are often 
modeled and interpreted in the context of self-organized crit¬ 
icality (SOC), where a nonlinear dynamical system reaches a 
stable critical point in which continuous energy input is re¬ 
leased intermittently through avalanches and in a scale-free 
way. SOC natural ly predicts power-laws in energy and WT 
distributions. See lAschwanden et all (l2014ll for a recent re¬ 
view on the many areas displaying SOC behavior. 

In astrophysics WTDs are studied for many d ifferent kinds 
of sources, such as out b ursting magnetars (iGogiis et al.l 
Il999t iG o&iis et aH l2000t iGavriil et akl l2004ft . flare stars 
(lArzner & Gudelll2004ft . and particularly the activity of the 
Sun throughout its cycle. WTDs of solar X-ray flares ex¬ 
hibit power-law tails with indices in the range 2.0-2.4 a cross 
several decades (iBoffetta et al.l[T9^ IWheatlandll2000h . de¬ 
pending on the class of flares and flux thresholds. Related 
bursty emission from the Sun such as coronal mass ejections 
(CMEs) are found to show very similar WTDs, whose in¬ 
dex ranges from ^1.9 to ^3.0 in low t o high-activity pe¬ 
riods of the solar cycl e (IW heatlandl l20()3h. L ikewise. WTDs 
of solar radio storms (lEastwood et al.ll2010ll . of solar ener¬ 
getic particle and o f solar electron events show very similar 
power-law indices (iLi et al.ll2014ft . Such power-law tailed 
WTDs are usually interpreted as the consequence of a time- 
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varying Poisson process produced by SOC systems, in which 
the energy input rate is intermitt ent and directly affects the 
degree of clusteriza tion of flares dAschwanden & McTiernanI 
l2010tlLr et alJl2014li . In this model, the bursty energy release 
is the result of avalanches produced in active regions where 
magnetic flux is twisted by the moving footpoints, leading 
to a series of independent magnetic reconnection events and 
consequent plasma acceleration. Alternatively to SOC, inter¬ 
pretations in the context of fully developed MHD turbulence 
have also been proposed to explain the bursty dynamics and 
the power-law WTD; the intermittent character is the result of 
a nonlinear dynamics which makes the convective motion of 
the fluid and magnetic field swing betwe en laminar and turbu¬ 
lent regimes repea tedly and chaotically (iBoffetta et al.ll 19^ 
iLepreti et alj|2004l) . 

The WTD between adjacent peaks in GRB 7 -ray prompt 
emission profiles was found to be describe d by a lognor¬ 
mal - which implies some degree of memory- (iLi & Fenimorel 
with an excess at long values due to QTs 

(iNakar & Piranl2002HOuilligan et al.l2002tlDrago & Pagliaral 

l2007h . However, when the peak detection efficiency is care¬ 
fully taken into account, it is found that the intrinsic WTD 
at short values is also compatible with an exponential, that 
is what i s expected for a constant Poi sson (thus memoryless) 
process (iBaldeschi & Guidorzill2015h . On the other side of 
the distribution long QTs could either mark the inner engine 
temporarily switching off, or result from modulatio n of the 
relativistic wind of shells (iRamirez-Ruiz et ^1200 Ih . or they 
could be due to a different physic al mechanism from that of 
short WTs (e.g.. lDrago et al.ir2008l) . 

In spite of the impressive data that are routinely being ac¬ 
quired in the Swift era, little progress has been reported on 
WTDs in GRBs. Recently, energy and WT distributions 
for GRB X-ray flares have been shown to have power-law 
tails very similar to what is observed for solar X-ray flares. 
In particular, the W TD has a power-law index of 1.8 ±0.2 
(IWang & Pail 1201^ . These results were interpreted as ev¬ 
idence for SOC possibly driven by magnetic reconnection 
episodes triggered in magnetized shells emitted by differen¬ 
tially rotating millisecond pulsars or, alternatively, by an h y- 
peraccreting disk around a black hole (iPopham et al.lll9^ . 

Yet, there are seyeral crucial issues which can be tackled 
with WTDs: is there additional eyidence for a link between 
prompt 7 -rays and late X-ray flares? To what extent do QTs 
differ from short WTs? Is it possible to proyide a common 
description of short WTs, QTs, X-ray flares? What about 
rest-frame properties? Is there eyidence for memory in GRB 
engines? What can be inferred on GRB engines through the 
WTD study? 

In this paper we address these issues through the analy¬ 
sis of the WTD of GRB prompt peaks for three independent 
data sets: Swift/BA^, CGRG/BATSE, and Fenni/GBM. For 
the Swift GRBs which haye also been promptly obseryed with 
XRT, we present a joint analysis of 7 -ray peaks and X-ray 
flares merged together. Section |2] describes the data sample 
selection and how we modeled the WTDs. Here we delib¬ 
erately did not consider the energy distribution of peaks and 
flares, because eyen though our peak search algorithm iden¬ 
tified moderately oyerlapped pulses, estimating their energy 
would require specific assumptions on their temporal struc¬ 
ture. We therefore decided to postpone it for future inyesti- 
gation. The results, their implications and interpretation are 
reported in Sections |3] and 0] respectiyely. Hereafter, uncer¬ 
tainties on best fit parameters are giyen at 90% confidence. 


unless stated otherwise. 

2. DATA ANALYSIS 

We se arched all long-durat ion 7 -ray light curyes with 
MEPS/Q (lGuidorzill2015l l^014ft . a peak search algorithm de¬ 
signed and calibrated to this goal. The adyantage of MEPSA 
compared with analo gous algorithms such as the one by 
ILi & Fenimorel (Il996h (LF) is twofold: 

• it has a lower false positiye rate. This is particularly true 
for the time interyals in which the signal drops to back¬ 
ground between two adjacent actiyity periods: in the 
best cases the LF false positiye rate is 3-5 x 10“^ bin“', 
while the MEPSA one is 1-2x10“^ bin“' jGuidorzil 

IMl; 

• it has a higher true positiye rate, especially at low signal 
to noise ratios ('^4-5). 

2 . 1 . Sample selection 

2.1.1. Swift/BAT data 

We started with the GRBs detected by SwiftfBAH in burst 
mode from January 2005 to September 2014, collecting 825 
GRBs. We extracted the mask-weighted light curyes in the 
15-150 keV energy band with a uniform bin time of 64 ms 
following the standard procedure recommended by the BAT 
tearr0 and applied MEPSA. We then imposed a minimum 
threshold of 5 cr significance, which ensures a yery low false 
positiye rate (< 10“^ bin“*; lGuidorzill2015h and selected the 
GRBs with at least two peaks. We then remoyed from our 
sample the short duration GRBs (both with and without ex¬ 
tended emission) by cross checking with the class ification pro- 
yided in the BAT catalog (ISakamoto et al.ll201 111 , as they will 
the subject of future inyestigation. Since this catalog does not 
include GRBs from 2010, for these GRBs we used the Tgo yal- 
ues as published in the BAT refined GCN circulars regularly 
published by the BAT team and set a conseryatiye minimum 
threshold of Tgo > 3 s. A couple of GRBs detected by BAT 
exhibited a yery long duration which could not be coyered en¬ 
tirely in burst mode: in one case we used the W/AD/Konus 
light curye for GRB 091024 (IVirgili et akllMl . while in the 
case of GRB 130 925A we us e d the peak times as they haye 
been obtained by lEyans et al.l (120141) from the corresponding 
Konus light curye. Finally, we ended up with a sample of 418 
long GRBs with at least two significant (> 5 (t) peaks each, 
totaling 2000 peaks and 1582 WTs. Hereafter, we refer to this 
sample as the BAT set. 

2.1.2. CGRO/BATSE data 

We took the concatenated 64-ms burst data distributed by 
the BATSF teamQ. For each curye we interpolated the back¬ 
ground by fitting with polynomials of up to fourth degree as 
suggested by the BATSF team ('e.g.. lGuidorzill2005h . Like in 
the BAT case, we applied MEPSA to an initial sample of 2024 
light curyes in the full passband. We applied the same se¬ 
lection on the peak significance and selected the long GRBs 
by requir ing Tgn > 2 s, where Tgo was taken from the GRB 
catalog (iPaciesas et al.llT99^ . We ended up with a sample 

^ littp://www.fe.infn.it/u/guidorzi/new_guidorzi_files/code.litml 
^ littp://swift.gsfc.nasa.gov/analysis/tlireads/bat_threads.html 
® ftp://cossc.gsfc.nasa.gov/compton/data/batse/ascii_data/64ms/ 

’ http://gammaray.msfc.nasa.gov/batse/grb/catalog/cun'ent/tables/duration_table.txt 


















































Waiting times in GRBs 


3 


of 1089 long GRBs with at least two 5 -(t significant peaks. 
Overall we collected 7649 peaks and 6560 WTs. Hereafter, 
this will be referred to as the BATSE sample. 

We also applied the same selection procedure to the light 
curves corresponding to the sum of the two softest energy 
channels (1 and 2 ) and to the sum of the two hardest channels 
(3 and 4), respectively within the 25-110 keV and >110 keV 
bands. We collected 1065 and 922 GRBs, with 5156 and 4912 
WTs, respectively. These two groups will be hereafter re¬ 
ferred to as BATSE 12 and BATSE34 sets, respectively. 

2.1.3. Fermi/GBM rfafa 

We selected 586 long GRBs detected with Eermi/GBM 
jMeegan et alj|200^ from July 2008 to December 2013. We 
extracted the light curves of the two brightest GBM units in 
the energy band 8-1000 keV with 64 ms resolution and sub¬ 
tracted the background through interpolation with a polyno¬ 
mial of up to third degree. We selected the long GRBs by 
imposing Tgo > 2 s, where 7go was taken from the official cat- 
alog0 We restricted to time intervals whose median values 
range from -30 to 300 s with reference to the trigger time. 
This corresponds to the time interval covere d by the time 
tagged event (TTE) data type in trigger mode (iPaciesas et alJ 
I 2 OI 2 I: iGruber et alj|2014b . Before -30 s and after 600 s the 
time resolution is that of CTIME data, 0.256 s. In most cases 
we did not consider intervals t > 300 s, because interpolation- 
estimated background often becomes critical an d the required 
effort for a proper estimate is beyond our scope (iGruber et alJ 
I 2 OI U lEitzpatrick et al.l l2012^ . We did not consider GRBs 
showing prolonged activity beyond this time interval. We then 
applied the same selection criteria as for the previous sets. 
Einally, by visual inspection we re moved phosphoresc ence 
spikes due to high-energy particles (iMeegan et al.ll200'^ . by 
comparing the same profiles in different GBM units. We 
ended up with a final sample of 2383 peaks out of 544 GRBs 
with at least two significant peaks. The total number of WTs 
is 1839. 

2.1.4. Swift/XRTX-rayflares 

We considered the catalog of 498 X -ray flare candidates 
detect ed with Swift/XRT obtained by ISwenson & Roming 
(l2014l) . This was extracted from 680 XRT light curves from 
January 2005 to December 2012 with a method based on 
the identification of breakpoints in the residuals of the fitted 
piecewise power-law light curves: these points mark sudden 
changes in the mean value due to unfitted features. The op¬ 
timal set of breakpoints was then found by minimizing the 
residual sum of squares against piecewise constant functions. 
To counter the effect of overfitting with unnecessary break- 
poin ts, they made use of t he Bayesian Information C riterion 
(see ISwenson et af] 1201 31 : ISwenson & Romin3l20I4l for fur¬ 
ther details). In this catalog each candidate is assigned a con¬ 
fidence value. We conservatively selected the subsample with 
a minimum confidence of 90%, ending up with 205 X-ray 
flare candidates. 

We separately merged each X-ray flare catalog with the 
BAT one by joining, for each GRB, the sequence of 7 -ray 
peak times and flare peak times into a unique sequence of tem¬ 
poral peaks. In doing this, every peak which was seen in both 
instruments was tagged as a 7 -ray peak and not considered 
any more as an X-ray flare. Analogously to the requirements 


for the previous sets, we selected those GRBs with at least 
two (either 7 or X-ray) peaks, so as to have at least one WT. 
We ended up with a sample of 1098 (954 7 - and 144 X-ray) 
peaks in 244 GRBs (01/2005 - 12/2012). We hereafter refer 
to this joint set as BAT-X sample. 

Einally, we selected the subsample with known redshift, so 
as to derive the WTD in the source rest-frame. This was done 
by simply correcting for cosmic dilation and thus dividing the 
observed WTs by the corresponding (1 -l-z). Unlike the width 
of a given pulse, which is affected not only by cosmic dilation 
but also by the energy passband shift, for their nature WTs 
are affected by the latter to a much lesser extent. We found 
359 WTs in 94 GRBs with known redshift. The subset with 
known redshift will be referred to as BAT-Xz. 

As an independent check, we in paralle l consi dered 
the X-ray flare catalo gs of IChincarini et al.l (1201 Oh and 
iBernardini et al.l (1201 Ih . which respectively include 113 
early-time (f < 10^ s) X-ray flares from April 2005 to March 
2008 and 36 late-time (f > 10^ s) flares from April 2005 to 
December 2009. However, due to lower statistics, we here¬ 
after focus on the BAT-X sample. 

2.2. Waiting time distribution modeling 

In physics a Poisson process is usually assumed to be char¬ 
acterized by a constant expected rate. The WTD of this pro¬ 
cess is exponential with e-folding r = 1 /A, 

P(Af) = = Ae-^^', (1) 

T 

where A is the constant mean rate and t is the mean WT. A 
time-varying Poisson process is characterized by a variable 
mean rate A(f): the process is locally Poisson, but the expected 
rate changes with time as described by A(f). According to 
this definition, the resulting process is the combination of two 
different proc esses at play and is o ften referred to as a “Cox 
process” (e.g.. lCox & Ishamlll980h : 

(a) : at a given time t events are generated according to a Pois¬ 

son process with rate A = A(f) and, as such, are statisti¬ 
cally independent; 

(b) : the expected rate A is itself a function of time, which 

can vary either randomly or deterministically as time 
passes. 

To derive the corresponding WTD, one may approximate 
A(f) as a piecewise constant function in a number of ad¬ 
jacent time intervals f,- (/=!,...,«) and treat it as a se- 
quence of several Poisson process es with rate A, . Poliowing 
Aschwanden & McTiernan| (l2010ll and references therein, the 
resulting WTD is 

P(Af) ^<^(A,)A, e-^- ^', (2) 

i 

where 

</.(A,) = , (3) 

is proportional to the expected number of WTs in interval f, 
where A = A, . In the continuous limit, Eq. (|2| becomes 

P(Af) = - 

/g A(f) dt 


http://heasarc.gsfc.nasa.gov AV3Browse/fermi/fermigbrst.html 


(4) 
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where T is the total duration. When A(f) is either unknown or 
hard to treat, it is possible to define /(A) such that f(X)dX = 
dt/T, that is the fraction of time during which the expected 
rate lies within the range [A, A + ^f A]. Equation (|4|i becomes, 


PijXt) = 


/;°°/(A)A^e-^^^^A 

/;“A/(A)t/A 


(5) 


We adopted the model for /(A) provided bv iLi et aT1(l2014h 
in their eq. (5), which has been proposed to fit the WTD 
obtained for solar X-ray flares and solar energetic particle 
events, 

/(A) = A A““ exp(-^ A), ( 6 ) 


with a and j3 free parameters and A is a normalization term 
(0 < Of < 2). This model generalizes several ot her models 
which had been put forward in the sam e context (IWheadandl 
l2000l:[Aschwanden & McTiernanll2010t) . The mean rate A is 

n+oo 

A = / Xf{X)dX=A/3°^-^r(2-a), (7) 

Jo 

where r() is the gamma function. From Equations (|5]l6]) the 
corresponding WTD is 


and it is normalized like a probability density function (pdf), 
i.e. /q P(At)d{At) = 1. There are only two free parameters, 
a, which determines the level of clusterization, and the char¬ 
acteristic WT P at which the WTD breaks; at At 2> /3, Eq. ® 
becomes a power-law with an index of (3 - a). 

Equation (| 6 ]l naturally gives rise to clusterization, i.e. time 
intervals characterized by an intense activity with a high rate 
of peaks (high A), interspersed with quiescent periods, during 
which the rate drops significantly (low A). The larger a, the 
shallower the power-law r egime at long WTs, and the more 
clustered the time profil e (lAschwanden & McTiernanI 120101: 
lAschwanden et al.1l2014l) . The details of how clustered the 
time profile looks like, in particular how the shot rate varies 
with time, are directly described by Eq. (HJ. At a given aver¬ 
age rate A, by increasing a the variance of A increases corre¬ 
spondingly, i.e. the shot rate varies more wildly. This means 
deviating more and more from the constant-rate case, thus en¬ 
hancing the clustering character. Figure [T] illustrates the dif¬ 
ference between a time-varying process like (a= 1, /3 = 0) in 
Eq. (| 6 ]l and a constant one sharing the same mean rate over a 
100-s time window. The temporal sequence of events for the 
former is evidently more clustered than that of the latter and, 
in spite of the typical fluctuations of a Poisson point process, 
tracks the behavior of A(f). It is worth nothing that, in gen¬ 
eral, in a Poisson process individual events are independent 
of each other and, as such, have no memory of the events that 
occurred earlier, regardless of whether the expected rate A is 
constant or time-dependent. The difference instead lies in the 
observed average rates as a function of time, so not on (a) but 
on (b); depending on whether A(f) varies either in a determin¬ 
istic way, or randomly with/without memory, the average rate 
inherits the corresponding degree of correlation. 

The WTDs we wanted to model are characterized by rare 
long WTs, where the count statistics is so low that one cannot 
use a simple minimization to fit the expected distribution 
of Eq. ([ 8 ]l to the counts collected in each bin. On the other 
side, merging the bins so as to have enough counts loses infor¬ 
mation and resolution. We therefore devised a log-likelihood 



Time [s] 

Fig. 1.— Example of time-varying Poisson mocess with variable rate A(0 
(dashed) assuming a = 1 and /3 = 0 in Eq. ^ and a constant one (solid) 
with the same mean values. Squares and vertical bars in the top region mai'k 
the corresponding event times that were generated as a consequence. The 
pronounced clusterization of the variable case over the constant one is clear. 



Fig. 2.— WTDs of 7 -ray pulses of the BAT (crosses), BATSE (circles), and 
GBM (squares) samples together with their corresponding best fit models. 
Error bars are the (normalized) square root of counts and are just indicative. 


based on Po isson statistics, which is essentially the C statistic 
(ICashlllQ^ and holds exactly even in the low count regime. 
We used it in the context of a Bayesian Markov Chain Monte 
Carlo approach. The details are reported in AppendixlAl 

3. RESULTS 

Table [1] reports the best fit parameters for all of the WTDs 
we considered. In all cases the model of Eq. (| 8 ]l provides an 
acceptable description. The lowest confidence level is that 
of BATSE (3.0%), still comparable with nominal 5% usu¬ 
ally adopted as a threshold. The BATSE sample is the largest 
(6560 WTs), so the high statistical sensitivity is likely to en¬ 
hance small deviations from the model. 

Figure displays the WTDs for the 7 -ray peak samples 
only. Apart from the GBM, whose power-law index is sig¬ 
nificantly steeper, the BAT and BATSE samples are fit with 
comparable indices, 2.067 og and 1.76 ±0.04, respectively. 
This is remarkable, given the different kind of detectors, en¬ 
ergy passbands, diffe rent GRB p opulations each instrument 
is mostly sensitive to (lBandll20^ . The soft and hard BATSE 
samples have the same index, showing no significant depen¬ 
dence on the energy channel. We investigated the reasons 
for the steeper WTD of the GBM set as follows; the dearth 
of long WTs is likely due to the short er scan ned time inter¬ 
vals, mostly from -30 to 300 s (Sect. 12.1.31 1. We therefore 
truncated the light curves of the Swift/BAH set and revised 
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TABLE 1 

Best fit parameters of the model in Eq. ((8]i obtained for 
DIFFERENT WTDS. SIZE IS THE NUMBER OF WTS. 


Sample 

Size 

a. 

V 

PL index 

^L 




(S) 

(=3-0:) 

(%) 

BAT 

BATSE 

1582 

6560 

0 94 +o.u'i 

1.24 ±0.04 

6 53'**''^^ 

2 06+® *® 

1.76 ±0.04 

26.4 

3.0 

BATSE 12 
BATSE34 

GBM 

BATtranc 

BAT-X 

BAT-Xz 

5156 

4912 

1839 

1445 

854 

359 

1.19 ±0.05 
1.18±0.05 
0.64)|';j« 

^•a‘+_o.o7 

1 AC-fO.lO 
-Q.ll 

'j 't'ttO.SS 

^■^^-0.16 
0./0_( 14 

(. 00+1.63 

6 33+*'^'^ 

20 

1 26+®''^^ 

1.81 ±0.05 

1.82±0.05 

2 36+®®'^ 

2 22^:1§ 
^•^^-0.15 

^•ot»_0,06 

- 0.10 

7.5 

76.6 

36.3 

5.2 

5.4 

18.5 


the WT selection accordingly. The results are reported in 
Table [T] as BATtrunc set, which includes 1445 WTs. Com¬ 
pared with the original BAT set, the WTD of the truncated 
data becomes steeper, from 2 . 06 ^j]o 9 to 2 . 22 ^q[ 5 , i.e. com¬ 
patible with the GBM value within uncertainties. Hence we 
interpret the slightly steeper value of the GBM set as the result 
of shorter time profiles which disfavor long WTs. 

Figure [3 displays the BAT-X set (squares) together with 
the corresponding best fit model. The power-law index is 
1.66;!^ Qg, i.e. compatible with BATSE sets within uncertain¬ 
ties. What is more, merging X-ray flares did not change the 
stochastic nature exhibited by the WTD, but extended its dy¬ 
namical range by at least one order of magnitude with WTs 
as long as 10 ^ s. A common stochastic model is found to well 
describe the WTD observed across more than five orders of 
magnitude. 

A similar result is obtained when one restricts to the 
known-redshift sample BAT-Xz in the GRB rest-frame (cir¬ 
cles in Fig. O: here the power-law index is 1.55 ;!;q [q, i.e. 
somewhat shallower. The rest-frame characteristic time is 
significantly shorter because of cosmic dilation, 1.3 s instead 
of the observer-frame values of 6-7 s. 

We also searched for possible correlations between WTs 
and peak intensities and between WTs and peak fluences of 
adjacent pulses, but found none. Finally, we repeated the 
analysis for various subsets of GRBs, by requiring a mini¬ 
mum number of of peaks per GRB and found no significant 
difference. 

3.1. y-ray peaks vs. X-ray flares 

Figure |4] shows the distribution of the number of 7 -ray 
peaks per GRB for two different classes of GRBs, depend¬ 
ing on whether their subsequent X-ray emission contains X- 
ray flares. Surprisingly, it is found that almost all GRBs 
(23/25) with A > 8 7 -ray peaks have no X-ray flares, al¬ 
though the two groups have comparable size, 131 and 163 
GRBs with and without flares, respectively. The two distribu¬ 
tions are unlikely to share a common population of events; a 
Kolmogorov-Smirnov (KS) test yields a mere 0.21% proba¬ 
bility, i.e. they are different with 3.1 ct (Gaussian) confidence. 
We visually inspected each of these 7 - and X-ray light curves 
and found only one case of a flareless GRB, whose X-ray 
light curve exhibited some low-level flaring activity which 
did not pass the 90%- confidence threshold in the flare sample 
selection (Sect. 12.1.41) . Therefore, GRBs with many pulses 
are unlikely to exhibit flares in the following declining X-ray 
emission. 

4. DISCUSSION 


Our results may be summarized in four fundamental as¬ 
pects: 

1. 7 -ray peaks and X-ray flares are compatible with being 
different aspects of the same stochastic process, which 
goes on after the end of the GRB itself and spans more 
than five orders of magnitude in time; 

2. short (interpulse) and long (quiescent) WTs between 7 - 
ray peaks are different realizations of the same stochas¬ 
tic process, the latters being only less frequent than the 
formers; hence a GRB with QTs and another without 
are by no means more different from each other than 
any other kind of GRBs are; 

3. GRBs with several (> 8 ) 7 -ray pulses are unlikely to 
exhibit X-ray flares after the end of the prompt emis¬ 
sion; 

4. 7 -ray peaks and X-ray flares tend to cluster in much the 
same way that solar flares, energetic particle events, and 
CMEs do, even though the processes may be different. 

The lognormal natu re of the WTD originally claimed 
(iFi & Fenimorel^~996^ has recently been shown to be possi¬ 
bly an artifact of the peak detection algorithms in the short 
WT end (< 1 s), w here peaks significantly over lap and can 
hardly be separated (iBaldeschi & Guidorzill2015h . We found 
that the long value (> few s) tail needs no more to be de¬ 
scribed as the sum of a lognormal tail and a power-law ex¬ 
cess due to the presen ce of QTs, that were interpreted as 
a different component (Nakar & Piranl 120021 lOuilligan et al.l 
l 2002 t IPragQ & PagliaralToOTh . ThiTapp^ent diversity also 
concerns the so-called precurso rs (lFazzatil2005HBurlon et al.l 
I 2 OO 8 I l200^ ICharisi et al.ll2014^ . which are nothing but emis¬ 
sion periods that are less intense than the following activity 
from which they are separated by a quiescent time. Our re¬ 
sults (1. and 2.) show that all kinds of WTs, including precur¬ 
sors, can be described within a common stochastic process, 
and this holds all the way up to late X-ray flares, thus point¬ 
ing towards a common mechanism, which keeps on working 
during and after the end of the prompt 7 -ray emission, before 
the afterglow emission due to the interaction with the external 
medium takes over. 

Another question concerns the break at low values in the 
WTD modeled in terms of the characteristic WT /3: is it an 
intrinsic property or is it entirely due to the low efficiency at 
short values of the peak detection algorithm? The capabil¬ 
ity of separating overlapping structures depends on a number 
of variables, such as the ratio between WT and peak widths, 
on intensities, and on temporal structures. While the drop at 
Af < 0.5 s is certainly due to the algorithm efficiency, the 
break itself modeled with /? is more complex: /3 is shorter 
at harder energies (Table [T]). A give n pulse has a narrowe r 
temporal structure at harder energies (iFenimore et al.l[r995h . 
whereas in t he softest energy ch annels there is a slow-varying 
component (IVetere et al.l 120061) . The presence of such soft 
component might hinder the peak identification in some case, 
so we examined the light curves in the harder channels. Vi¬ 
sual inspection suggests that the paucity of subsecond WTs 
with respect to the power-law extrapolation is real and is un¬ 
likely to be a mere artifact of the peak identification process. 
In addition, minimum pulse widths o bserved in GRB pro¬ 
files typically are i n the range 01-1 s (Fenimore et akllTO^ 
iNorris et all 119^ iMarguttiet all 1201 ll) . The MEPSA effi¬ 
ciency is above 10% for such pulse widths, for WTs > 0.5 s. 
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Fig. 3.— WTDs of 7 -ray pulses and X-ray flares of BAT-X (squares) and of BAT-Xz (circles) sets together with their con'esponding best fit models. Eri'or bars 
are the (normalized) square root of counts and are just indicative. 



Fig. 4.— Distribution of number of 7 -ray peaks per GRB for two distinct 
subsets of the Swift BAT-X sample, depending on the presence/lack of flares 
in the following X-ray emission. Almost all (23/25) GRBs with at least 8 
7 -ray peaks have no X-ray flares. A KS test yields a common population 
probability of 0 . 21 %. 

and for measured signal to noise ratios > 5 (lGuidorzill2015l) . 
It is therefore unlikely that the algorithm efficiency is entirely 
responsible for the observed exponential cutoff observed in 
the low end of the WTDs. 

4.1. A simple toy model 

We devised a very simple toy model to explore more 
in detail how a time-dependent Poisson process like the 
one of Eq. (| 6 ]l could be obtained in a GRB engine. For 
the sake of clarity, suppose each pulse marks the accre¬ 
tion of a single fragment of an hyperaccreting disk. Actu¬ 
ally, the idea behind this model is more general and only 
deals with the sequence of bursty emission episodes and 
their probability of occurring within a given time; how¬ 
ever, hereafter we refer to the model of an hyperaccret¬ 


ing disk being fragmented as the source of the stochastic 
proce ss which is re s ponsible for the prompt 7 -r av emis- 


i ss which IS re s ponsible tor the prompt 7 -r av emis- 
(lWooslevlll993l: iMacFadven & Woosle\^ll999h as well 
for the subsequent X-ray flare activ ity ( King_etaflj200^ 
Perna et al.l 1200(4 iKumar et alJ 120081 : iCannizzo & Gehrelsl 
2009tlGeng et aklbool) . We briefly summarize the basic in- 


sion 

as 


gredients of the model, which are then thoroughly described 
in the remaining part of this Section: 


• a number of fragments are independently accreted with 
the same, constant, probability per unit time; 


• the number of available fragments is obviously decreas¬ 
ing with time; this naturally leads to a time-dependent 
Poisson process whose mean accretion rate decreases 
with time; 


• at the beginning, if the mean rate is too high (A > l/P), 
accretion becomes inefficient and is suppressed by a 
factor of exp(-A/3); 

• for some (~ 30%) GRBs the reservoir of fragments is 
split into two separate groups sharing the same individ¬ 
ual accretion probability per unit time, but with the sec¬ 
ond group becoming available only at later times (e.g., 
the late group could be identified with the outer part of 
the accretion disk). 

Fet us assume that the disk or the inner part of it has been 
split into a number of fragments, each of which has the same 
given probability of accreting per unit time, independently of 
the others. The probability for a given fragment to survive up 
to a given time t is proportional to exp(-f/r), where r is the 
mean accretion time for each fragment. The total expected 
rate scales as the number of fragments that are still available, 
A = |iV(f)| = (No/t) expi-t/r) =N{t)/T. The analogous /(A) 
to Eq. © is found as 

/(A)iiA oc |^f^| =r ^ , (9) 
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At[s] 

Fig. 5.— Series of exponential WTDs (thin solid) of a sequence of con¬ 
stant Poisson processes with progressively increasing e-folding. At interme¬ 
diate values the mean WTD (thick solid) scales as a power-law with index 2, 
shown for comparison (dashed). 

which corresponds to the a = 1 case in Eq. (| 6 ]) at A < 1//3. 
Rather than a continuously varying, A(f) of this model is de¬ 
scribed more exactly by a piecewise Poisson process like the 
one of Eqs. dull, where A, = i/r is the expected rate when 
i fragments are left over. All terms have equal weights ^,’s, 
since each piecewise constant process contributes one WT. 
The resulting WTD is thus given by Eq. dll). 


P(Af) 


, No . 
No 

l=l 


( 10 ) 


which can also be expressed as, 

_ .r [Nox^^^^-(No+ 1 ] 

Not(I-x)^ 


( 11 ) 


where x = _ In Figure |5] an example of such WTD is 

shown, with No = 20 initial fragments, r = 1 s. As time goes 
by, A(f) decreases and the e-folding of the individual expo¬ 
nential WTDs (thin solid) increase correspondingly. As a re¬ 
sult, the total WTD (thick solid) show a power-law regime 
with index 2 at intermediate values of Af. At At < t /No 
the total WTD is dominated by the initial exponential with 
e-folding t/Nq, when fragments are all available. This 
agrees with the result that the intrinsic (i.e., corrected for 
the algorithm efficiency) WTD at short values is likely ex- 
ponential, that is, compatible with a constant Poisson process 
ilBaldeschi & Guidorzill201^ . 

Thus far, with reference to Eq. dH) our model implicitly as¬ 
sumed 13 = 0 (see Eq. |9]). However, in our attempt to repro¬ 
duce the observed WTD with the piecewise constant process 
of Eq. dTOl) failed to model the observed break at Af ~ 1//3. 
So we required that, when the expected rate becomes com¬ 
parable or higher than 1//3, the number of observed WTs 
is suppressed by a factor of exp(-A^) with respect to our 
model. This can be interpreted as if, when the number N 
of fragments that can be readily accreted is such that the ex¬ 
pected rate is A = N/t > 1//3, the overall process becomes 
inefficient and the rate is suppressed by exp(-A^). In other 
words, the number of WTs shorter than /3 is smaller than 
what is expected from Eq. dU. This introduces some degree 
of memory in the initial stages of the accretion process: as 
long as the number of fragments ready to be accreted is too 
high (t/N < /3) some of them are temporarily halted from 



Fig. 6. — WTD for the Swift BAT-X sample (squares) compai'ed with a 
simulated sample of 903 WTs derived from a toy model (circles). The shaded 
interval is where the peak search algorithm efficiency drops. 


accreting by some mechanism connected with the accretion 
rate itself. For instance, this self-regul ating mechanism could 
be driven by the magnetic held (e.g., iProga & Zhan jl2006l : 

Uzdenskv & MacFadvenl2006l : lBernardini et alJ20131) 7 which 

is known to have a complex role in accretion processes of ut¬ 
terly different objects such as T Tauri stars (IStenhens et alJ 
12014) . However, we cannot provide a more specihc and 
physical justihcation for the exponential character of this 
self-quenching mechanism, which is therefore ad-hoc in its 
present formulation. 

We assumed the logarithmic average and l-cr scatter of the 
BAT-X WTD, 16.8 s and a multiplicative scatter of 7.1, to 
generate the values for r for each simulated burst. The num¬ 
ber of generated peaks in each simulated curve was taken 
from the observed distribution and was augmented by 20 % 
to ensure that the detected peaks were enough (since some are 
missed by the algorithm). The peak times for each simulated 
curve were randomly generated from an exponential distribu¬ 
tion with e-folding r, i.e. independently from each other. To 
mimic the drop in the peak detection efficiency at short WTs 
as well as the mechanism mentioned above about the suppres¬ 
sion at high rates, we overlooked each peak occurring within 
0.5 s of the previous one, while through a binomial we as¬ 
signed each At a probability exp(-/3/Af) of being observed, 
where /? was set to the htted value of the real WTD (Table [1]). 

To obtain a good match with the observed WTD over the 
same range we had to make a further assumption: we as¬ 
sumed that for a fraction of GRBs (~ 30%) randomly selected 
through a binomial, the disk is fragmented equally in two 
groups, the hrst of which is available for being accreted from 
the beginning (to = 0 ), while the second one becomes available 
from fo = 50 T on, where r is the common mean accretion time 
for each individual fragment from t = to- The reason behind 
this is the observation of two similar bunches of 7 -ray peaks 
interspersed with a long quiescent time (up to several tens of 
seconds) for a small fraction of GRBs. Physically, this could 
be the result of an outer part of the disk being accreted at later 
times with respect to the inner one or, more in general, de¬ 
layed additional energy reservoir becoming available for late 
internal dissipation, with minimum variability tim escale com¬ 
parab l e with that of the early prompt emission (iFan & Weil 
120051 iLazzati & Pernall20d^ iTroia et akllMT^ . While the 
choice of the fraction of such GRBs and the duration of the 
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Fig. 7.— Distribution of the ratio of adjacent WTs of the Swift BAT-X 
sample (solid) and of the toy-model sample (shaded). 

quiescence period are somewhat arbitrary, the good match be¬ 
tween simulated and real WTD does not depend crucially on 
them. Overall, the goal here is just to show the plausibility 
of the essential properties of this model, which can reproduce 
the observed properties in spite of the simple assumptions. We 
ended up with a set of 903 simulated WTs, whose distribution 
is compared with the real one in Fig.| 6 ] 

We further tested the consequences of this toy model by 
studying the distribution of the ratio between adjacent WTs. 
While WTDs describe how WTs distribute as a whole, losing 
information on their temporal sequence, the ratio distribution 
focuses on that. We therefore selected from the BAT-X as 
well as from the simulated sample the GRBs with > 3 peaks, 
so as to have at least two WTs, and derived the two distribu¬ 
tions shown in Fig. \J\ A KS test between the two sets yields 
43% probability that they were drawn from a common popu¬ 
lation. Logarithmic mean and dispersion for the real (simu¬ 
lated) data are p. = 0.14 and cr = 0.72 (p = 0.09 and a = 0.81). 
Similar results are obtained adopting other non-parametric 
tests, such as the Wilcoxon -Mann-Whitney or the m ore sen¬ 
sitive Epps-Singleton one (lEpps & Singletonlll986h . respec¬ 
tively yielding 45% and 9% probability. Interestingly, simply 
replacing Eq. (|9]l with a constant Poisson process and apply¬ 
ing the very same following steps, one ends up with a nar¬ 
rower and more zero-centered logarithmic ratio distribution, 
p = 0.013 and a = 0.32, which is rejected with high confi¬ 
dence (p-value < 2 X 10“'®) from a KS test. This means that 
for a constant process the ratio is, on average, closer to 1 , and 
is less scattered around it than real data. The compatibility of 
the ratio distribution predicted by the toy model with the real 
one proves that the temporal sequence of WTs is compatible 
with an evolving Poisson process and is incompatible with a 
constant one on long timescales. In particular. X-ray flares 
are nothing but some of the last fragments that are left over 
and which are accreted on long timescales, when the rate de¬ 
creases in a granular way, following the very same stochastic 
process ruling the accretion of the earliest ones. Flence, no 
correlation is to be expected between 7 -ray prompt emission 
duration (T gn) and X-ray flar es times, in agreement with ob- 
servations (ILiang et al.ll2006l) . Einally, the result of Sect. 13.11 
can be easily explained: GRBs with many 7 -ray peaks ac¬ 
crete fragments rapidly with relatively short r, so that at late 
time very few or none at all are left over for X-ray flares. 
The same result could be explained differently though: multi- 
peaked GRBs could have on average a brighter early X-ray 


afterglow continuum which overshines possible X-ray flares, 
which would then go unseen. 

Overall, we assumed a direct connection between emission 
and observed times of the peaks. Within the context of in¬ 
ternal shocks the observed time profile of both prompt 7 - 
rays and lat e X-ray flares is strictly connected to the emis - 
sion history (iKobavashi et al.ll997HMaxham & Zhangl2009h . 
Should this not be the case, little could be inferred about on 
emission times -and potentially the times of individual accre¬ 
tion episodes- from the study of the observed WTD. However, 
this connection becomes more complex due to the variety of 
Lorentz factors associated with the wind of shells colliding 
with each other. Even though the intrinsic duration of the 
GRB engine a ctivity may differ by a factor of a few from the 
observed one (iGao & Meszarosll2014l) . on average the tempo¬ 
ral sequence of mutual collisions between randomly-assigned 
Lorentz-factor shells should track the emission time history. 
The nature of a given WTD is not altered as a whole when 
one passes from the emission to the observed times: in fact, 
each shell has a Lorentz factor -which is in principle what can 
make the observed WTs very different from the emitted ones- 
that is independent of the WTs preceding and following that 
shell. This statistical independence ensures that the observed 
WTD keeps memory on the emission time distribution. Only 
at late times, when the average Lorentz factor is expected to 
systematically decrease and the statistical independent char¬ 
acter likely begins to fail, long WTs are likely to be affected 
as a consequence. 

4.2. Solar activity: analogies and differences 

It is remarkable and intriguing that WTDs of both solar 
eruptive events (X-ray flares, radio storms, high-energy par¬ 
ticle events, CMEs) and of GRBs can be modeled with the 
same kind of time-dependent Poisson process. The power- 
law characterization of the WTD heavy tail must not be over¬ 
interpreted from a mathematical viewpoint, since power-laws 
are, in general, what one ends up with when dealing with the 
sum of independent heavy-tailed variables. It works much in 
the same way that a normal distribution is the final outcome 
of the sum of independent finite-variance variables. In ad¬ 
dition, claiming that data are power-law distributed is con¬ 
trived whene ver the explored range does not cover at least 
two decades (IStumpf & Ported l2012h . In this sense, invok¬ 
ing a SOC-driven mechanism for GRBs purely based on 
the power-law character of the WTD, and possibly of the 
energy distribution too, as suggested for X-ray flares from 
GRBs (IW ang & Daill2013h or from other black hole systems 
(IWang et al.ll20l4 . is a stretched interpretation of the data, as 
we explain below. 

The same or very similar power-law indices imply that both 
processes have very similar degrees of clusterization, with 
analogous swings between intense and low-activity periods, 
apart from temporal rescaling (seconds for GRBs, hours for 
the Sun). However, one has to be careful extending this anal¬ 
ogy to a common physical mechanism. Overall, there is a fun¬ 
damental difference in terms of dynamical systems between 
GRB inner engines during core collapse and the Sun: for 
the latter, the regions where eruptive phenomena take place 
continuously receive energy, which is then released through 
avalanches, thus making the SOC interpretation plausible (al¬ 
though alternatives based on MHD turbulence seem equally 
compatible with observations). Instead, GRB engines are sys¬ 
tems which start with a very-far-from-equilibrium configu¬ 
ration, evolve very fast using up all of the available energy. 
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which -no matter how much- is limited. A GRB inner engine 
cannot return to its original configuration, it goes through an 
obviously irreversible evolution, whereas this is not the case 
for the solar active regions over sufficiently long timescales. 
For this reason, one needs not invoke SOC mechanisms re¬ 
lated to accretion disks, in particular there is no need for a 
mechanism like the one pro posed to explain 1 / f fl uctuations 
in black hole power spectra (iMineshige et al.l[l994l) . 

Therefore, a simple time-varying Poisson process explains 
the secular evolution of the mean rate of bursts/flares as well 
as the stochastic independent character of individual energy 
release episodes. This model disregards the physical origin 
of fragmentation and how energy is distributed among differ¬ 
ent fragments: thus, in principle it is co mpatible with var ious 
physical drivers, such as gra vitational (|Perna et al.ll2006l) . or 
magnetorotational instability (iProga & Zhang^l2006^ . 

5. CONCLUSIONS 

In this paper we studied the waiting time distributions 
of GRB 7 -ray pulses in three catalogs, CGRG/BATSE, 
Fermi/GUM, and Swift/BAH. For the latter, for the first time 
we merged 7 -ray pulses and X-ray flares detected with 
Swift/XRT belonging to the same GRBs, and for a subsam¬ 
ple the same analysis was carried out in the source rest- 
frame. We found that all WTDs can be described in terms of 
a common time-varying Poisson process that rules different 
waiting time intervals, which thus far in the literature have 
been treated differently: specifically, we showed that short 
WTs (< 1 s), long quiescent times (> 10 s) all the way up 
to late time X-ray flares are the manifestation of a common 
stochastic process. GRB WTDs exhibit heavy tails which are 
modeled with power-laws over 4-5 decades in time with in¬ 
dices in the range 1 .7-2.1, depending on the relative weight 
of late time events, such as X-ray flares, in each GRB sam¬ 
ple. Because of the ubiquitous nature of power-laws (cen¬ 
tral limit theorem for heavy-tailed distributions), the charac¬ 
ter of WTDs must not be imbued with a mystical sense or 
overinterpreted as evidence for a universal process. In this 
sense, the similarity of the WTD power-law index with that 
of solar eruptive phenomena, such as flares and coronal mass 
ejections, proves nothing but a similar degree of clusteriza¬ 
tion in time. Nonetheless, it is remarkable that the WTD of 
7 -ray pulses and that of X-ray flares not only have compat¬ 
ible power-law indices but they join and extend the dynami¬ 
cal range for a common sample of GRBs. All this points to 
a common stochastic process ruling both phenomena. The 


unification under a common process of all different kinds of 
waiting times in GRBs (short interpulse times, long quiescent 
times, time intervals following precursors, time intervals be¬ 
tween the end of the 7 -ray prompt emission and subsequent 
X-ray flares) is a new result. 

Another noteworthy result is that GRBs with many (> 8 ) 7 - 
ray pulses are unlikely (3.1a confidence in Gaussian units) to 
exhibit X-ray flares in their subsequent early X-ray emission. 
This result is naturally explained in the context of the time- 
varying Poisson process: many pulses observed in the prompt 
of a given GRB are indicative of a relatively short mean ac¬ 
cretion time for a single disk fragment. Consequently, most 
of the available fragments are consumed during the prompt 
phase, with very few or none at all left over for the subsequent 
phase. 

In the light of the irreversible evolution of GRB inner en¬ 
gines, the interpretation of a time-varying Poisson process ap¬ 
pears to be simple and reasonable: the secular evolution of 
the expected rate of events is naturally linked to the energy 
reservoir being gradually used up, whereas the stochastically 
independent accretion of individual fragments is explained by 
the Poissonian character of the process. 

Although self-organized criticality models naturally pre¬ 
dict power-law tailed distributions of waiting times and en¬ 
ergy, drawing upon this kind of dynamics for GRBs might be 
premature. Other equally plausible alternatives, such as fully 
developed MHD turbulence, can explain the same properties, 
as it was suggested for the solar case. Possible evidence for 
turbulence in GRBs has also been suggested from the analy¬ 
sis of power density spectra (iBeloborodov et alJll998[ 120001 : 
iGuidorzi et al.ll2012l : lDichiara et al.ll2013h . Yet. we find that a 

simple time-varying Poisson process such as that of a system 
gradually using up all the available pieces already provides a 
remarkably accurate description. 

The energy distribution, which was beyond the scope of 
this paper, will help to further constrain the stochastic process 
and possibly clarify whether more complex dynamical mod¬ 
els, such as SOC or MHD turbulence, are to be considered. 

We are grateful to the anonymous referee for a constructive 
and insightful review. PRIN MIUR project on “Gamma Ray 
Bursts: from progenitors to physics of the prompt emission 
process”, P. I. F. Frontera (Prot. 2009 ERC3HT) is acknowl¬ 
edged. 


APPENDIX 


LOG-LIKELIHOOD TO FIT THE DISTRIBUTIONS 

Let the WTD consist of N logarithmically spaced bins, each collecting C, counts. Let Af, 1 and Af, 2 be the lower and upper 
bounds of the /-th bin (i = 1,...,N). Integrating Eq. ® within this time interval yields the corresponding expected counts, 
Ei(a, P), where we explicitly meant that it depends on the model parameters: 


/-Ar;; 


Ei(a,P) = Ctc 


P(Af) t/(Af) = C,o 


ol-a 




(/3 + Af,-,i)“-2-(/3+Af,-,2)' 


ol-2 


(Al) 


where Ctot = Ym=\ /3 is a function of both model parameters (Sect. lA2l l. The probability of C, counts is ruled by the Poisson 

distribution where £, is the expected value. 


P(C,|a,/3) = e-E- 



(A2) 
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where the dependence on model parameters is implicit through The total probability is thus 


P(C\a,p) 


He- 


Q! ’ 


(A3) 


where C is the set of observed counts per bin {C, } (i = 1,... ,N). The corresponding negative log-likelihood is therefore 


L{a,l3) 


N N 

-^log(P(C,|a,/3)) = ^(£, + log(C,!)-C, log£,) . 

?=1 ?=1 


(A4) 


We determine the best fit model parameters and their uncertainties in the Bayesian context. The posterior probability density 
function of the parameters for a given observed distribution C, is (Bayes theorem) 


P(a,l3\C) 


P(C\a,f3)Pia,f3) 

PiC) 


(A5) 


where the first term in the numerator of the right-hand side of eq. (IA51 l is the likelihood function of Eq. (IA31 l. Pia, (3) is the prior 
distribution of the model parameters, and the denominator is the normalization term. We assumed a uniform p rior distribution, 
since no a priori informati on is available on the model parameters. The mode of the posterior probability of Eq. (IA5I) is therefore 
found by minimizing Eq. (IA4t . 

We estimate the posterior density of the model parameters through a Markov Chain Monte Carlo (M CMC) algorithm such 
as the random-walk Metropolis-Hastings in the implementation of the package MHadaptiv43 (v.1.1-8). We initially 
approximate the posterior using a bivariate normal distribution centred on the mode and with covariance matrix obtained by 
minimizing Eq. (IA4b . Eor each WTD we generate 5.1 x 10"^ sets of simulated model parameters and retain one every 5 MCMC 
iterations after excluding the first 1000. The remaining 10"^ sets of parameters are therefore used to approximate the posterior 
density. Einally, once the best fit model parameters are determined, the bivariate posterior distribution of (q;,/3) is sampled via 
MCMC simulations, which yield expected value and 90% confidence intervals for each of them. 

As a matter of fact, since the bins in the low end of distribution are strongly affected by the poor efficiency of MEPSA, these 
are to be ignored. In practice, one has to replace in Eq. (lAlb Ctot with where k\ and k 2 are the first and last bins to 

be considered. In addition, the same £, in Eq. (lAlb has to be further divided by a renormalizing factor so that it becomes. 


Ei(a,l3) 


(/3 + Af,-,i)"-^-(/3 + Af,-,2)"-" 

(/3+A4..i)“-2-(/3+Af,„2)“-2 


(A6) 


Eor the WTDs discussed in the present paper we considered Af > 0.5 s (Af > 0.2 s) in the observer (source) rest frame. 

To assess the goodness of the fit for a given WTD, we use each set of simulated values for (a, 13) to generate as many synthetic 
WTDs from the the posterior predictive distribution. Hence for a given observed WTD, we directly calcu late 10"^ synthetic 
realizations of the same WTD. Eor each of these WTDs we calculate the negative log-likelihood with Eq. (IA4b and derive a 
corresponding distribution of values, against which the value obtained from the real WTD is checked. This comparison directly 
provides a confidence level of modelling the observed WTD in terms of the best fit model of Eq. ([8]l. 
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